## Callin Switzer
## 17 Nov 2017
## Multilevel model to visualize bees'
## behavior on the artificial pollen system
#install packages
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c("ggplot2", "reshape2", 'lme4', 'sjPlot', "multcomp", "plyr", "effects")
ipak(packages)
## Loading required package: ggplot2
## Loading required package: reshape2
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: sjPlot
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.10
## Current Matrix version is 1.2.11
## Please re-install 'TMB' from source or restore original 'Matrix' package
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
## Loading required package: plyr
## Loading required package: effects
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
## ggplot2 reshape2 lme4 sjPlot multcomp plyr effects
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
# set ggplot theme
theme_set(theme_classic() + theme(axis.text=element_text(colour="black")))
# define data and figure directories
dataDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehData"
figDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehFigs"
print(paste("last run ", Sys.time()))
## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'zone/tz/2017c.1.0/
## zoneinfo/America/Los_Angeles'
## [1] "last run 2017-12-08 22:58:40"
print(R.version)
## _
## platform x86_64-apple-darwin15.6.0
## arch x86_64
## os darwin15.6.0
## system x86_64, darwin15.6.0
## status
## major 3
## minor 4.2
## year 2017
## month 09
## day 28
## svn rev 73368
## language R
## version.string R version 3.4.2 (2017-09-28)
## nickname Short Summer
Read in data and double check it
sl <- read.csv(file.path(dataDir, '01_CombinedTrials_cleaned.csv'))
# make sure hive is a factor
sl$hive = as.factor(sl$hive)
head(sl)
## beeCol hive meanFreq IT_imputed index freq amp
## 1 blue 4 395.2778 3.61 1 400 0.06328
## 2 blue 4 395.2778 3.61 2 360 0.34299
## 3 blue 4 395.2778 3.61 3 430 0.44099
## 4 blue 4 395.2778 3.61 4 420 0.16841
## 5 blue 4 395.2778 3.61 5 420 0.15284
## 6 blue 4 395.2778 3.61 6 410 0.18302
## datetime rewNum rewTF lowRewAmp highrewAmp
## 1 2016_11_22__08_23_47_721 1 T 0 5
## 2 2016_11_22__08_23_49_677 2 T 0 5
## 3 2016_11_22__08_23_50_139 3 T 0 5
## 4 2016_11_22__08_24_10_938 4 T 0 5
## 5 2016_11_22__08_24_18_013 5 T 0 5
## 6 2016_11_22__08_24_18_464 6 T 0 5
## BeeNumCol
## 1 BeeBlue_22Nov2016_Hive4_initial
## 2 BeeBlue_22Nov2016_Hive4_initial
## 3 BeeBlue_22Nov2016_Hive4_initial
## 4 BeeBlue_22Nov2016_Hive4_initial
## 5 BeeBlue_22Nov2016_Hive4_initial
## 6 BeeBlue_22Nov2016_Hive4_initial
## accFile trialNum
## 1 2016_11_22__08_23_47_721_220_450_test.txt 1
## 2 2016_11_22__08_23_49_677_220_450_test.txt 1
## 3 2016_11_22__08_23_50_139_220_450_test.txt 1
## 4 2016_11_22__08_24_10_938_220_450_test.txt 1
## 5 2016_11_22__08_24_18_013_220_450_test.txt 1
## 6 2016_11_22__08_24_18_464_220_450_test.txt 1
## datetime_str lowFrq highFrq IT trt amp_acc
## 1 2016-11-22 08:23:47.721000 220 450 3.61 full 6.222222
## 2 2016-11-22 08:23:49.677000 220 450 3.61 full 33.725664
## 3 2016-11-22 08:23:50.139000 220 450 3.61 full 43.361849
## 4 2016-11-22 08:24:10.938000 220 450 3.61 full 16.559489
## 5 2016-11-22 08:24:18.013000 220 450 3.61 full 15.028515
## 6 2016-11-22 08:24:18.464000 220 450 3.61 full 17.996067
dim(sl) # should be 24303 rows
## [1] 24303 21
# hive 5 is most common
table(sl$hive, useNA = 'always')
##
## 3 4 5 <NA>
## 513 1047 22743 0
## make sure all bee colors are lowercase
sl$beeCol <- tolower(sl$beeCol)
# make sure there are values lower than 220 and higher than 450
# (the cutoff for buzzes used in the experiment)
hist(sl$freq, breaks = seq(215, 450, by = 10))
nrow(sl[sl$freq < 220 | sl$freq > 450,]) # should have 0 rows
## [1] 0
# look at treatments
xtabs(~sl$beeCol+ trt, data = sl )
## trt
## sl$beeCol full high low
## blue 36 0 0
## gold 54 0 0
## goldred 3 0 0
## green 13 0 0
## lime 28 0 0
## limeblue 530 0 0
## limegold 54 0 0
## limegreen 50 0 0
## limeorange 84 0 0
## limepink 51 0 0
## limepurple 109 0 243
## limepurpleyellow 58 0 0
## limered 50 0 3424
## limesilver 3 0 0
## limewhite 101 0 0
## limeyellow 52 0 0
## orange 85 416 0
## orangeblue 50 0 0
## orangegreen 50 0 0
## orangepink 50 0 0
## orangepurple 50 0 0
## purple 9 0 0
## redblue 673 0 0
## redgreen 530 0 0
## redpink 56 0 621
## redpurple 55 163 0
## silver 26 0 0
## white 283 0 0
## whiteblue 56 0 2730
## whitegold 48 0 1177
## whitegreen 39 0 0
## whiteorange 54 0 1068
## whitepink 135 1049 0
## whitepurple 55 0 0
## whitered 59 1285 0
## whiteyellow 157 0 0
## yellowblue 54 1533 0
## yellowgreen 532 0 0
## yelloworange 54 2059 0
## yellowpink 58 0 2251
## yellowpurple 54 0 1189
## yellowred 547 0 0
# find percentage reward by treatment
mean(grepl("[tT]", as.character(sl$rewTF))) # overall mean
## [1] 0.4155454
# percentage that were rewarded by treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, mean)
## full high low
## 1.0000000 0.3535742 0.2128631
# total number of trials for each treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, length)
## full high low
## 5095 6505 12703
# total number of rewards per treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, FUN = function(x) sum(x))
## full high low
## 5095 2300 2704
# total number of trials that were unrewarded per treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, FUN = function(x) sum(!x))
## full high low
## 0 4205 9999
Visualizations and modeling
# calculate trial averages and plot
sl$colNum = paste(sl$beeCol, sl$hive, sep = "_")
sl$trt <- as.character(sl$trt)
sl$trt[sl$trt == "full" & sl$trialNum >1 ] <- "full_2"
aggdata <- aggregate(sl$freq, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(aggdata)[colnames(aggdata) == "x"] = "freq"
aggdata_sd <- aggregate(sl$freq, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=sd, na.rm=TRUE)
colnames(aggdata_sd)[colnames(aggdata_sd) == "x"] = "freq_sd"
aggdata = merge(aggdata, aggdata_sd)
#aggdata$trt = sapply(1:nrow(aggdata), FUN = function(x){sl[sl$colNum == aggdata$colNum[x] & sl$trialNum == aggdata$trialNum[x], "trt"][1]})
aggdata = aggdata[order(aggdata$colNum, aggdata$trialNum, decreasing = FALSE), ]
agg_sm = aggdata[aggdata$trialNum <= 2, ]
rownames(agg_sm) = 1:nrow(agg_sm)
agg_sm
## colNum trialNum trt freq freq_sd
## 1 blue_4 1 full 395.2778 39.38717
## 2 gold_3 1 full 348.3333 37.15089
## 3 goldred_4 1 full 406.6667 32.14550
## 4 green_4 1 full 318.8889 57.32461
## 5 green_4 2 full_2 285.0000 59.16080
## 6 lime_5 1 full 374.6429 36.36055
## 7 limeblue_5 1 full 302.4000 34.73280
## 8 limeblue_5 2 full_2 299.8148 39.02059
## 9 limegold_5 1 full 344.6296 44.07168
## 10 limegreen_5 1 full 370.8000 29.40498
## 11 limeorange_5 1 full 285.7576 35.26953
## 12 limeorange_5 2 full_2 344.5098 26.70683
## 13 limepink_5 1 full 308.4314 43.23760
## 14 limepurple_5 1 full 352.2642 38.66232
## 15 limepurple_5 2 low 344.2798 32.55948
## 16 limepurpleyellow_5 1 full 336.7241 50.76219
## 17 limered_5 1 full 354.6000 37.91559
## 18 limered_5 2 low 365.3381 36.20344
## 19 limesilver_5 1 full 306.6667 11.54701
## 20 limewhite_5 1 full 292.0000 44.90352
## 21 limewhite_5 2 full_2 327.0588 35.51305
## 22 limeyellow_5 1 full 339.4231 39.37818
## 23 orange_3 1 full 388.5294 38.22837
## 24 orange_5 1 full 345.0980 33.60789
## 25 orangeblue_5 1 full 301.4000 45.17585
## 26 orangegreen_5 1 full 296.2000 40.45103
## 27 orangepink_5 1 full 286.6000 27.15113
## 28 orangepurple_5 1 full 329.0000 25.65469
## 29 purple_3 1 full 382.2222 21.08185
## 30 redblue_4 1 full 342.2449 40.53113
## 31 redblue_4 2 full_2 335.6667 51.89189
## 32 redgreen_5 1 full 334.7059 34.19666
## 33 redgreen_5 2 full_2 314.5283 31.53545
## 34 redpink_5 1 full 288.3929 41.06654
## 35 redpink_5 2 low 275.3448 34.24450
## 36 redpurple_5 1 full 311.4545 31.64726
## 37 redpurple_5 2 high 336.4444 40.75995
## 38 silver_5 1 full 305.3846 38.39070
## 39 white_4 1 full 343.1250 47.49720
## 40 white_4 2 full_2 342.2388 38.95689
## 41 whiteblue_5 1 full 318.2143 27.17667
## 42 whiteblue_5 2 low 351.6242 29.13402
## 43 whitegold_5 1 full 315.2083 34.20772
## 44 whitegold_5 2 low 307.2727 45.76009
## 45 whitegreen_4 1 full 300.7692 33.43372
## 46 whiteorange_5 1 full 324.4444 51.71280
## 47 whiteorange_5 2 low 355.2308 44.74234
## 48 whitepink_5 1 full 336.8750 42.68285
## 49 whitepink_5 2 high 323.9161 28.55770
## 50 whitepurple_5 1 full 328.1818 43.50781
## 51 whitered_5 1 full 288.8136 36.05956
## 52 whitered_5 2 high 304.6667 42.10035
## 53 whiteyellow_5 1 full 303.4545 39.07344
## 54 whiteyellow_5 2 full_2 312.5532 36.56252
## 55 yellowblue_5 1 full 313.7037 44.86021
## 56 yellowblue_5 2 high 327.9688 49.76099
## 57 yellowgreen_5 1 full 298.0392 40.64577
## 58 yellowgreen_5 2 full_2 300.6000 27.58364
## 59 yelloworange_5 1 full 268.1481 36.18996
## 60 yelloworange_5 2 high 313.9205 39.62646
## 61 yellowpink_5 1 full 323.6207 39.98676
## 62 yellowpink_5 2 low 335.7042 41.00507
## 63 yellowpurple_5 1 full 354.4444 36.37678
## 64 yellowpurple_5 2 low 373.0058 34.17237
## 65 yellowred_5 1 full 320.1754 41.85396
## 66 yellowred_5 2 full_2 333.0357 29.96481
agg_sm[duplicated(data.frame(agg_sm$colNum, agg_sm$trialNum)), ]
## [1] colNum trialNum trt freq freq_sd
## <0 rows> (or 0-length row.names)
ggplot(agg_sm, aes(x = trt, y = freq, fill = trialNum > 1)) +
geom_boxplot(alpha = 0.2) +
geom_point(aes(color = trialNum>1)) +
geom_line(aes(group = colNum))
diffdf <- sapply(unique(agg_sm$colNum), FUN = function(x){
tmp = agg_sm[agg_sm$colNum == x, ]
if(nrow(tmp) <= 1)
diff = NA
else
diff = tmp$freq[tmp$trialNum == 2] - tmp$freq[tmp$trialNum == 1]
return(diff)
})
trtDF = sapply(unique(agg_sm$colNum), FUN = function(x){
tmp = agg_sm[agg_sm$colNum == x, ]
ttrs = paste(tmp$trt[tmp$trialNum == 1], tmp$trt[tmp$trialNum == 2], sep = "_")
return(ttrs)
})
buzzdiffs = data.frame(trtDF, diffdf)
ggplot(buzzdiffs, aes(x = trtDF, y= diffdf)) +
geom_boxplot() +
geom_point()
## Warning: Removed 20 rows containing non-finite values (stat_boxplot).
## Warning: Removed 20 rows containing missing values (geom_point).
agg2 = aggregate(sl$freq, by=list(colNum = sl$colNum, fullTrt = sl$trt == "full", trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(agg2)[colnames(agg2) == "x"] = "freq"
agg2$trt = as.character(agg2$trt)
diffdf <- t(as.data.frame(t(sapply(unique(agg2$colNum), FUN = function(x){
tmp = agg2[agg2$colNum == x, ]
if(nrow(tmp) <= 1)
return(NA)
if (length(unique(tmp$trt)) > 2){
tmp = tmp[tmp$trt != "full_2", ]
}
diff = tmp$freq[!tmp$fullTrt] - tmp$freq[tmp$fullTrt]
return(diff)
}))))
length(diffdf)
## [1] 43
trtDF = sapply(unique(agg2$colNum), FUN = function(x){
tmp = agg2[agg2$colNum == x, ]
if (length(unique(tmp$trt)) > 2){
tmp = tmp[tmp$trt != "full_2", ]
}
ttrs = paste(tmp$trt[tmp$fullTrt], tmp$trt[!tmp$fullTrt], sep = "_")
return(ttrs)
})
length(trtDF)
## [1] 43
buzzdiffs = data.frame(trtDF, diffdf)
tapply(buzzdiffs$diffdf, INDEX = buzzdiffs$trtDF, mean)
## full_ full_full_2 full_high full_low
## NA 4.851689 13.209001 14.307123
ggplot(droplevels(buzzdiffs[buzzdiffs$trtDF != "full_", ]), aes(x = trtDF, y= as.numeric(diffdf))) +
geom_boxplot() +
geom_point() +
labs(y = "diff in frequency -- averages")
m1 = lm(as.numeric(diffdf) ~ trtDF - 1, data = droplevels(buzzdiffs[buzzdiffs$trtDF != "full_", ]))
summary(m1)
##
## Call:
## lm(formula = as.numeric(diffdf) ~ trtDF - 1, data = droplevels(buzzdiffs[buzzdiffs$trtDF !=
## "full_", ]))
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.741 -17.992 -3.556 15.767 53.901
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## trtDFfull_full_2 4.852 7.488 0.648 0.524
## trtDFfull_high 13.209 9.667 1.366 0.186
## trtDFfull_low 14.307 8.372 1.709 0.102
##
## Residual standard error: 23.68 on 21 degrees of freedom
## Multiple R-squared: 0.1987, Adjusted R-squared: 0.08422
## F-statistic: 1.736 on 3 and 21 DF, p-value: 0.1904
sl$trt = relevel(factor(sl$trt), ref = "full")
sl$trt2 = mapvalues(sl$trt, from = c("full", "full_2", "high", "low"),
to = c("full", "full", "high", "low"))
# summary for paper
# fit a varying slope and intercept for colNum (bee ID), and allow the slope of the trialNum
# variable to vary by colNum (beeID)
m2 = lmer(freq ~ trt2 + hive + trialNum + IT_imputed + (1+trialNum|colNum), data = sl, REML = FALSE)
# this model estimates a global intercept
# random effect intercept for colNum (beeID)
# a single global estimate for trialNum
# the effect of trialNum within each level of colNum
# the correlation between intercept of trialNum across levels of colNum
m2_1 = lmer(freq ~ trt2* IT_imputed + hive + trialNum + (1+trialNum|colNum), data = sl, REML = FALSE)
# compare models with BIC
BIC(m2, m2_1) # m2 (no interaction) is better
## df BIC
## m2 11 246086.5
## m2_1 13 246090.2
anova(m2, m2_1) # Note: this disagrees with BIC
## Data: sl
## Models:
## m2: freq ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum |
## m2: colNum)
## m2_1: freq ~ trt2 * IT_imputed + hive + trialNum + (1 + trialNum |
## m2_1: colNum)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2 11 245997 246087 -122988 245975
## m2_1 13 245985 246090 -122979 245959 16.547 2 0.0002552 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 = update(m2, .~. - hive)
BIC(m2, m3) # m3 is better (without hive)
## df BIC
## m2 11 246086.5
## m3 9 246073.9
anova(m2, m3) # Note: disagrees with BIC
## Data: sl
## Models:
## m3: freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum | colNum)
## m2: freq ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum |
## m2: colNum)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 9 246001 246074 -122991 245983
## m2 11 245997 246087 -122988 245975 7.5417 2 0.02303 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit best model with REML = TRUE
m3 <- update(m3, .~., REML = TRUE)
summary(m3) # summary for paper
## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum | colNum)
## Data: sl
##
## REML criterion at convergence: 245962
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0242 -0.5161 0.1663 0.6809 3.9135
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## colNum (Intercept) 1072 32.74
## trialNum 112 10.58 -0.73
## Residual 1439 37.93
## Number of obs: 24303, groups: colNum, 43
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 479.138 44.361 10.801
## trt2high 11.161 2.294 4.865
## trt2low 17.718 1.916 9.248
## trialNum 1.448 2.133 0.679
## IT_imputed -37.291 10.773 -3.462
##
## Correlation of Fixed Effects:
## (Intr) trt2hg trt2lw trilNm
## trt2high 0.069
## trt2low -0.005 0.009
## trialNum 0.007 -0.042 -0.037
## IT_imputed -0.993 -0.072 0.000 -0.093
plot(m3)
qqnorm(ranef(m3)$colNum[[1]])
qqline(ranef(m3)$colNum[[1]])
sjp.lmer(m3) # plot random effects to find any outliers
## Plotting random effects...
sjp.lmer(m3, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## Computing p-values via Wald-statistics approximation (treating t as Wald z).
# post-hoc tests -- pvals for paper
summary(glht(m3, linfct = mcp(trt2 = "Tukey")), test = adjusted("none"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum |
## colNum), data = sl, REML = TRUE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## high - full == 0 11.161 2.294 4.865 1.14e-06 ***
## low - full == 0 17.718 1.916 9.248 < 2e-16 ***
## low - high == 0 6.557 2.975 2.204 0.0275 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
Generate CI’s
# set number of bootstrap samples
nbootSims = 1000
ITmean = mean(tapply(sl$IT_imputed, INDEX = sl$beeCol, FUN = function(x) x[1] ))
# don't need hive, because that's not in the model we chose (above)
pframe <- data.frame(trt2 = levels(droplevels(sl$trt2)), IT_imputed = ITmean, colNum = 99999, trialNum = 2)
pframe$freq <- 0
pp <- predict(m3, newdata = pframe, re.form=NA, type = 'response') # re.form sets all random effects to 0
### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m3, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 1000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-bb2_se[1,]
pframe$bhi<-bb2_se[2,]
pframe$predMean <- pp
pframe <- pframe[, c('trt2', "blo", "bhi", "predMean")]
pframe
## trt2 blo bhi predMean
## 1 full 322.2094 336.6227 329.4150
## 2 high 332.2079 349.1315 340.5761
## 3 low 339.6179 354.8610 347.1330
Make plots for paper
#plot 95% confidence intervals
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
# rename labels for plot
pframe$trt3 = mapvalues(pframe$trt2, from = c("full", "high", "low"),
to = c("Full range\n(220 - 450 Hz)", "High range\n(340 - 390 Hz)", "Low range\n(220 - 330 Hz)"))
pframe
## trt2 blo bhi predMean trt3
## 1 full 322.2094 336.6227 329.4150 Full range\n(220 - 450 Hz)
## 2 high 332.2079 349.1315 340.5761 High range\n(340 - 390 Hz)
## 3 low 339.6179 354.8610 347.1330 Low range\n(220 - 330 Hz)
g0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
geom_point()+
labs(y = "Sonication frequency (Hz)", x = "Frequency range for reward") +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 360, label=c("a", "b", "c"),
color="black")
g0
ggsave(plot = g0, filename = file.path(figDir, "SonicationFreqPredsAndCI_unadjusted.pdf"), width = 5, height = 4)
g0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
geom_point()+
labs(y = "Sonication frequency (Hz)", x = "Frequency range for reward") +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 360, label=c("a", "b", "b"),
color="black")
g0
ggsave(plot = g0, filename = file.path(figDir, "SonicationFreqPredsAndCI_adjustedPvals.pdf"), width = 5, height = 4)
# use the "effects" library to get CI's
# this basically gives the same result as bootstrapping
# it's a good double check
ee <- as.data.frame(Effect(c("trt2"),
fixed.predictors =list(given.values = c(trialNum = 2,
IT_imputed = ITmean)),
m3) )
# compare two methods -- very similar
ee
## trt2 fit se lower upper
## 1 full 329.4150 3.798940 321.9689 336.8612
## 2 high 340.5761 4.280727 332.1856 348.9666
## 3 low 347.1330 4.081495 339.1330 355.1330
pframe
## trt2 blo bhi predMean trt3
## 1 full 322.2094 336.6227 329.4150 Full range\n(220 - 450 Hz)
## 2 high 332.2079 349.1315 340.5761 High range\n(340 - 390 Hz)
## 3 low 339.6179 354.8610 347.1330 Low range\n(220 - 330 Hz)
g1 <- ggplot(ee, aes(x=trt2, y=fit))+
geom_point()+
labs(y = "Sonication Frequency (Hz)", x = "Treatment") +
geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1)+
theme_classic() +
annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 360, label=c("a", "b", "c"),
color="black")
g1
Analysis for amplitude
# refref: here try log amp_acc and REML = FALSE with BIC
# interaction model
# note log-transformation to make model fit assumptions better
maa0 = lmer(log(amp_acc) ~ trt2* IT_imputed + hive + trialNum + (1+trialNum|colNum), data = sl, REML = FALSE)
# main effect model
maa1 = lmer(log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1+trialNum|colNum), data = sl, REML = FALSE)
BIC(maa0, maa1) # use no interaction (keep maa1)
## df BIC
## maa0 13 44377.70
## maa1 11 44360.43
anova(maa0, maa1) # note that this agrees with likelihood ratio test
## Data: sl
## Models:
## maa1: log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum |
## maa1: colNum)
## maa0: log(amp_acc) ~ trt2 * IT_imputed + hive + trialNum + (1 + trialNum |
## maa0: colNum)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## maa1 11 44271 44360 -22125 44249
## maa0 13 44272 44378 -22123 44246 2.9209 2 0.2321
maa2 = update(maa1, .~. - trt2)
BIC(maa1, maa2) # keep treatment (maa1)
## df BIC
## maa1 11 44360.43
## maa2 9 44375.90
anova(maa1, maa2) # agrees with LRT; p-val for paper, if needed
## Data: sl
## Models:
## maa2: log(amp_acc) ~ hive + trialNum + IT_imputed + (1 + trialNum |
## maa2: colNum)
## maa1: log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum |
## maa1: colNum)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## maa2 9 44303 44376 -22142 44285
## maa1 11 44271 44360 -22125 44249 35.665 2 1.801e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
maa3 <- update(maa1, .~. - hive)
BIC(maa1, maa3) # remove hive (maa3)
## df BIC
## maa1 11 44360.43
## maa3 9 44345.85
anova(maa1, maa3) # agrees with LRT
## Data: sl
## Models:
## maa3: log(amp_acc) ~ trt2 + trialNum + IT_imputed + (1 + trialNum |
## maa3: colNum)
## maa1: log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum |
## maa1: colNum)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## maa3 9 44273 44346 -22128 44255
## maa1 11 44271 44360 -22125 44249 5.6143 2 0.06038 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit with REML = TRUE for paper
maa3 <- update(maa3, .~., REML = TRUE)
summary(maa3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ trt2 + trialNum + IT_imputed + (1 + trialNum |
## colNum)
## Data: sl
##
## REML criterion at convergence: 44279.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8377 -0.5961 0.1156 0.7367 2.9119
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## colNum (Intercept) 0.07079 0.26607
## trialNum 0.00203 0.04506 -0.57
## Residual 0.35872 0.59893
## Number of obs: 24303, groups: colNum, 43
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.398060 0.455225 5.268
## trt2high 0.194526 0.034446 5.647
## trt2low 0.067759 0.029326 2.311
## trialNum 0.005513 0.010107 0.545
## IT_imputed 0.324122 0.110268 2.939
##
## Correlation of Fixed Effects:
## (Intr) trt2hg trt2lw trilNm
## trt2high 0.107
## trt2low -0.009 0.029
## trialNum 0.036 -0.089 -0.083
## IT_imputed -0.995 -0.114 0.000 -0.087
# diagnostics
plot(maa3)
qqnorm(ranef(maa1)$colNum[[1]])
qqline(ranef(maa1)$colNum[[1]])
sjp.lmer(maa3) # plot random effects to find any outliers
## Plotting random effects...
sjp.lmer(maa3, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## Computing p-values via Wald-statistics approximation (treating t as Wald z).
# post-hoc tests -- pvals for paper
summary(glht(maa3, linfct = mcp(trt2 = "Tukey")), test = adjusted("none"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(amp_acc) ~ trt2 + trialNum + IT_imputed +
## (1 + trialNum | colNum), data = sl, REML = TRUE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## high - full == 0 0.19453 0.03445 5.647 1.63e-08 ***
## low - full == 0 0.06776 0.02933 2.311 0.02086 *
## low - high == 0 -0.12677 0.04459 -2.843 0.00447 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
# post-hoc tests with bonf adjustment
summary(glht(maa3, linfct = mcp(trt2 = "Tukey")), test = adjusted("bonf"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(amp_acc) ~ trt2 + trialNum + IT_imputed +
## (1 + trialNum | colNum), data = sl, REML = TRUE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## high - full == 0 0.19453 0.03445 5.647 4.89e-08 ***
## low - full == 0 0.06776 0.02933 2.311 0.0626 .
## low - high == 0 -0.12677 0.04459 -2.843 0.0134 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)
Generate CI’s for amplitude
# set number of bootstrap samples
nbootSims2 = 1000
# don't need to include hive, but will set trialNum to 2
pframe <- data.frame(trt2 = levels(droplevels(sl$trt2)),
IT_imputed = ITmean, colNum = 99999, trialNum = 2)
pframe$amp_acc <- 0
# exponentiate to get back to original scale
pp <- exp(predict(maa3, newdata = pframe, re.form=NA, type = 'response')) # re.form sets all random effects to 0
### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(maa3, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims2)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 1000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
#exponentiate to get back to orignal scale
pframe$blo<-exp(bb2_se[1,])
pframe$bhi<-exp(bb2_se[2,])
pframe$predMean <- pp
pframe <- pframe[, c('trt2', "blo", "bhi", "predMean")]
pframe
## trt2 blo bhi predMean
## 1 full 39.05177 45.01982 41.91285
## 2 high 46.44376 56.04234 50.91299
## 3 low 41.17250 48.98075 44.85125
Make plots for amplitude for paper
#plot 95% confidence intervals
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
pframe$trt3 = mapvalues(pframe$trt2, from = c("full", "high", "low"),
to = c("Full range\n(220 - 450 Hz)", "High range\n(340 - 390 Hz)", "Low range\n(220 - 330 Hz)"))
pframe
## trt2 blo bhi predMean trt3
## 1 full 39.05177 45.01982 41.91285 Full range\n(220 - 450 Hz)
## 2 high 46.44376 56.04234 50.91299 High range\n(340 - 390 Hz)
## 3 low 41.17250 48.98075 44.85125 Low range\n(220 - 330 Hz)
ga0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
geom_point()+
labs( y = expression ("Sonication amplitude "(m~s^{-2})), x = "Frequency range for reward") +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = 'none') +
theme_classic() +
annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 56, label=c("a", "b", "c"),
color="black")
ga0
ggsave(plot = ga0, filename = file.path(figDir, "SonicationAmpPredsAndCI_unadjusted.pdf"), width = 5, height = 4)
# make plot using bonferroni adjustments
ga0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
geom_point()+
labs( y = expression ("Sonication amplitude "(m~s^{-2})), x = "Frequency range for reward") +
geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = 'none') +
theme_classic() +
annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 56, label=c("a", "b", "a"),
color="black")
ga0
ggsave(plot = ga0, filename = file.path(figDir, "SonicationAmpPredsAndCI_BonfAdjusted.pdf"), width = 5, height = 4)
Visualize aggregated data for amplitude
aggdata <- aggregate(sl$amp_acc, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(aggdata)[colnames(aggdata) == "x"] = "amp"
aggdata_sd <- aggregate(sl$amp_acc, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=sd, na.rm=TRUE)
colnames(aggdata_sd)[colnames(aggdata_sd) == "x"] = "amp_sd"
aggdata = merge(aggdata, aggdata_sd)
aggdata = aggdata[order(aggdata$colNum, aggdata$trialNum, decreasing = FALSE), ]
agg_sm = aggdata[aggdata$trialNum <= 2, ]
rownames(agg_sm) = 1:nrow(agg_sm)
agg_sm
## colNum trialNum trt amp amp_sd
## 1 blue_4 1 full 40.03641 27.005004
## 2 gold_3 1 full 32.58012 13.395404
## 3 goldred_4 1 full 27.49328 8.575236
## 4 green_4 1 full 44.14520 45.654277
## 5 green_4 2 full_2 19.19985 21.543257
## 6 lime_5 1 full 72.53866 33.940021
## 7 limeblue_5 1 full 49.20279 24.837998
## 8 limeblue_5 2 full_2 59.04592 36.444346
## 9 limegold_5 1 full 34.95839 26.390545
## 10 limegreen_5 1 full 52.56061 32.949195
## 11 limeorange_5 1 full 43.28056 22.966005
## 12 limeorange_5 2 full_2 50.75113 25.907474
## 13 limepink_5 1 full 49.57493 19.711454
## 14 limepurple_5 1 full 43.98249 22.702402
## 15 limepurple_5 2 low 39.95722 19.359774
## 16 limepurpleyellow_5 1 full 34.54496 17.151702
## 17 limered_5 1 full 42.27701 20.963674
## 18 limered_5 2 low 40.63711 20.895896
## 19 limesilver_5 1 full 26.24910 7.336471
## 20 limewhite_5 1 full 37.65233 22.862123
## 21 limewhite_5 2 full_2 40.91118 19.936507
## 22 limeyellow_5 1 full 43.51579 23.156058
## 23 orange_3 1 full 32.25152 18.932874
## 24 orange_5 1 full 55.08235 20.197053
## 25 orangeblue_5 1 full 51.97426 23.301810
## 26 orangegreen_5 1 full 58.59866 33.410963
## 27 orangepink_5 1 full 50.98635 28.370137
## 28 orangepurple_5 1 full 34.53312 16.739467
## 29 purple_3 1 full 53.61466 21.935890
## 30 redblue_4 1 full 28.88506 11.312274
## 31 redblue_4 2 full_2 48.61796 25.277613
## 32 redgreen_5 1 full 78.72927 46.273733
## 33 redgreen_5 2 full_2 70.38830 47.080737
## 34 redpink_5 1 full 74.20294 37.536798
## 35 redpink_5 2 low 56.85081 31.951114
## 36 redpurple_5 1 full 58.06343 40.818616
## 37 redpurple_5 2 high 53.13637 26.813039
## 38 silver_5 1 full 46.34517 30.280907
## 39 white_4 1 full 47.76100 21.162920
## 40 white_4 2 full_2 46.45015 22.564532
## 41 whiteblue_5 1 full 66.11499 33.532056
## 42 whiteblue_5 2 low 67.00735 29.330027
## 43 whitegold_5 1 full 43.10847 26.148044
## 44 whitegold_5 2 low 45.20369 27.241996
## 45 whitegreen_4 1 full 38.87762 18.184245
## 46 whiteorange_5 1 full 61.63802 39.835205
## 47 whiteorange_5 2 low 45.85399 25.838373
## 48 whitepink_5 1 full 63.53760 35.030181
## 49 whitepink_5 2 high 80.19250 47.237615
## 50 whitepurple_5 1 full 47.77247 35.753931
## 51 whitered_5 1 full 48.37233 25.860700
## 52 whitered_5 2 high 74.79446 45.764818
## 53 whiteyellow_5 1 full 65.50425 32.808944
## 54 whiteyellow_5 2 full_2 56.00923 33.372943
## 55 yellowblue_5 1 full 52.45834 30.688932
## 56 yellowblue_5 2 high 79.36156 41.580725
## 57 yellowgreen_5 1 full 61.58710 33.835105
## 58 yellowgreen_5 2 full_2 79.18968 27.697874
## 59 yelloworange_5 1 full 73.32616 36.109449
## 60 yelloworange_5 2 high 54.48161 24.329362
## 61 yellowpink_5 1 full 52.74199 24.869065
## 62 yellowpink_5 2 low 57.76271 29.939244
## 63 yellowpurple_5 1 full 62.77375 29.437162
## 64 yellowpurple_5 2 low 68.98565 35.865129
## 65 yellowred_5 1 full 62.77866 39.753892
## 66 yellowred_5 2 full_2 47.97187 28.376661
agg_sm[duplicated(data.frame(agg_sm$colNum, agg_sm$trialNum)), ]
## [1] colNum trialNum trt amp amp_sd
## <0 rows> (or 0-length row.names)
ggplot(agg_sm, aes(x = trt, y = amp, fill = trialNum > 1)) +
geom_boxplot(alpha = 0.2) +
geom_point(aes(color = trialNum>1)) +
geom_line(aes(group = colNum))
diffdf <- sapply(unique(agg_sm$colNum), FUN = function(x){
tmp = agg_sm[agg_sm$colNum == x, ]
if(nrow(tmp) <= 1)
diff = NA
else
diff = tmp$amp[tmp$trialNum == 2] - tmp$amp[tmp$trialNum == 1]
return(diff)
})
trtDF = sapply(unique(agg_sm$colNum), FUN = function(x){
tmp = agg_sm[agg_sm$colNum == x, ]
ttrs = paste(tmp$trt[tmp$trialNum == 1], tmp$trt[tmp$trialNum == 2], sep = "_")
return(ttrs)
})
buzzdiffs = data.frame(trtDF, diffdf)
ggplot(buzzdiffs, aes(x = trtDF, y= diffdf)) +
geom_boxplot() +
geom_point()+
labs(y = "amplitude difference m/s/s")
## Warning: Removed 20 rows containing non-finite values (stat_boxplot).
## Warning: Removed 20 rows containing missing values (geom_point).
agg2 = aggregate(sl$amp_acc, by=list(colNum = sl$colNum, fullTrt = sl$trt == "full", trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(agg2)[colnames(agg2) == "x"] = "amp"
agg2$trt = as.character(agg2$trt)
diffdf <- t(as.data.frame(t(sapply(unique(agg2$colNum), FUN = function(x){
tmp = agg2[agg2$colNum == x, ]
if(nrow(tmp) <= 1)
return(NA)
if (length(unique(tmp$trt)) > 2){
tmp = tmp[tmp$trt != "full_2", ]
}
diff = tmp$amp[!tmp$fullTrt] - tmp$amp[tmp$fullTrt]
return(diff)
}))))
length(diffdf)
## [1] 43
trtDF = sapply(unique(agg2$colNum), FUN = function(x){
tmp = agg2[agg2$colNum == x, ]
if (length(unique(tmp$trt)) > 2){
tmp = tmp[tmp$trt != "full_2", ]
}
ttrs = paste(tmp$trt[tmp$fullTrt], tmp$trt[!tmp$fullTrt], sep = "_")
return(ttrs)
})
length(trtDF)
## [1] 43
buzzdiffs = data.frame(trtDF, diffdf)
tapply(buzzdiffs$diffdf, INDEX = buzzdiffs$trtDF, mean)
## full_ full_full_2 full_high full_low
## NA 2.871714 9.194234 3.249068
ggplot(droplevels(buzzdiffs[buzzdiffs$trtDF != "full_", ]), aes(x = trtDF, y= as.numeric(diffdf))) +
geom_boxplot() +
geom_point()+
labs(y = "amplitude difference m/s/s")
Plot amplitude vs. frequency for initial trial, accounting for size
ggplot(sl[sl$trialNum == 1, ], aes(x = freq, y = amp_acc)) +
geom_point(position = position_jitter(width = 3), alpha = 0.2) +
theme_classic() +
geom_smooth(aes(group = trt2)) +
labs( y = expression ("Sonication amplitude "(m~s^{-2})), x = "Frequency (Hz)")
## `geom_smooth()` using method = 'gam'
centerFreq = scale(sl[sl$trialNum == 1, "freq"])
m1 <- lmer(log(amp_acc)~centerFreq + I(centerFreq^2) + I(centerFreq^3) + IT_imputed + (1|beeCol), data = sl[sl$trialNum == 1, ])
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ centerFreq + I(centerFreq^2) + I(centerFreq^3) +
## IT_imputed + (1 | beeCol)
## Data: sl[sl$trialNum == 1, ]
##
## REML criterion at convergence: 3526.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5361 -0.6344 0.0693 0.6857 2.4239
##
## Random effects:
## Groups Name Variance Std.Dev.
## beeCol (Intercept) 0.08395 0.2897
## Residual 0.32877 0.5734
## Number of obs: 1971, groups: beeCol, 42
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.00589 0.57016 3.518
## centerFreq 0.37124 0.03004 12.357
## I(centerFreq^2) -0.08590 0.01227 -7.003
## I(centerFreq^3) -0.03835 0.01055 -3.636
## IT_imputed 0.43530 0.13776 3.160
##
## Correlation of Fixed Effects:
## (Intr) cntrFr I(F^2) I(F^3)
## centerFreq -0.047
## I(cntrFr^2) -0.093 0.081
## I(cntrFr^3) -0.018 -0.844 0.065
## IT_imputed -0.996 0.043 0.071 0.019
plot(m1)
# predict -- using mean IT span
ppdf <- data.frame(centerFreq = sort(unique(centerFreq)), IT_imputed = ITmean, beeCol = 99999, acc_amp = 0)
# exponentiate to get back to original scale
ppdf$predAmp = exp(predict(m1, newdata = ppdf, type = "response", re.form = NA))
cf_Unscaled = ppdf$centerFreq * attr(centerFreq, 'scaled:scale') + attr(centerFreq, 'scaled:center')
ggplot(ppdf, aes(x = cf_Unscaled, y = predAmp))+
geom_line() +
labs(x = "Sonication Frequency (Hz)", y = expression ("Sonication amplitude "(m~s^{-2})) ) +
geom_point(data = sl[sl$trialNum == 1, ],
aes(x = freq, y = amp_acc),
alpha = 0.2, position = position_jitter(width =2), pch =16, size = .5)
ggsave(filename = file.path(figDir, "FreqVsAmp_1stTrial_rawData.pdf"), width = 4, height = 3)
g22 <- ggplot(ppdf, aes(x = cf_Unscaled, y = predAmp))+
geom_line() +
labs(x = "Sonication Frequency (Hz)", y = expression ("Predicted sonication amplitude "(m~s^{-2})) )
g22
ggsave(filename = file.path(figDir, "FreqVsAmp_1stTrial.pdf"), width = 4, height = 3)
sessionInfo()
## R version 3.4.2 (2017-09-28)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] effects_4.0-0 carData_3.0-0 plyr_1.8.4 multcomp_1.4-8
## [5] TH.data_1.0-8 MASS_7.3-47 survival_2.41-3 mvtnorm_1.0-6
## [9] sjPlot_2.4.0 lme4_1.1-14 Matrix_1.2-11 reshape2_1.4.2
## [13] ggplot2_2.2.1
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-131 pbkrtest_0.4-7 RColorBrewer_1.1-2
## [4] rprojroot_1.2 tools_3.4.2 TMB_1.7.11
## [7] backports_1.1.1 R6_2.2.2 sjlabelled_1.0.5
## [10] DT_0.2 lazyeval_0.2.1 mgcv_1.8-20
## [13] colorspace_1.3-2 nnet_7.3-12 tidyselect_0.2.3
## [16] mnormt_1.5-5 compiler_3.4.2 sandwich_2.4-0
## [19] labeling_0.3 scales_0.5.0 lmtest_0.9-35
## [22] psych_1.7.8 blme_1.0-4 stringr_1.2.0
## [25] digest_0.6.12 foreign_0.8-69 minqa_1.2.4
## [28] rmarkdown_1.8 stringdist_0.9.4.6 pkgconfig_2.0.1
## [31] htmltools_0.3.6 pwr_1.2-1 htmlwidgets_0.9
## [34] rlang_0.1.4 shiny_1.0.5 bindr_0.1
## [37] zoo_1.8-0 dplyr_0.7.4 magrittr_1.5
## [40] modeltools_0.2-21 bayesplot_1.4.0 Rcpp_0.12.13
## [43] munsell_0.4.3 abind_1.4-5 prediction_0.2.0
## [46] stringi_1.1.6 yaml_2.1.14 merTools_0.3.0
## [49] snakecase_0.5.1 grid_3.4.2 parallel_3.4.2
## [52] sjmisc_2.6.2 forcats_0.2.0 lattice_0.20-35
## [55] ggeffects_0.2.2 haven_1.1.0 splines_3.4.2
## [58] sjstats_0.12.0 knitr_1.17 codetools_0.2-15
## [61] stats4_3.4.2 glue_1.2.0 evaluate_0.10.1
## [64] modelr_0.1.1 nloptr_1.0.4 httpuv_1.3.5
## [67] gtable_0.2.0 purrr_0.2.4 tidyr_0.7.2
## [70] assertthat_0.2.0 mime_0.5 coin_1.2-1
## [73] xtable_1.8-2 broom_0.4.2 survey_3.32-1
## [76] coda_0.19-1 tibble_1.3.4 arm_1.9-3
## [79] glmmTMB_0.1.4 bindrcpp_0.2